#!/usr/bin/python

"""

  #--------------------------------------------------------------------------
  # pyCoord
  # pyCoord.py
  #
  # (c) 2010 Andy Turner
  # (c) 2005 Jonathan Stott
  #
  # This is a python translation of Jonathan Stott's PHPcoord, which was
  # made available under the terms of GPLv2.
  # 
  # Please refer to http://www.jstott.me.uk/phpcoord/ for history etc of the
  # original code.
  #--------------------------------------------------------------------------


  # ================================================================== LatLng

"""
  
from math import cos, sin, acos, asin, floor, atan, tan
  
class LatLng:
    
    lat = 0
    lng = 0
    
    """
     * Create a new LatLng object from the given latitude and longitude
     *
     * @param lat latitude
     * @param lng longitude
    """
    def __init__(self, lat, lng):
        self.lat = lat
        self.lng = lng

    """
     * Return a string representation of this LatLng object
     *
     * @return a string representation of this LatLng object
    """
    def __str__(self):
        return "(%f,%f)" % (self.lat, self.lng)
    
    """
     * Calculate the surface distance between this LatLng object and the one
     * passed in as a parameter.
     *
     * @param to a LatLng object to measure the surface distance to
     * @return the surface distance
    """
    def distance(self, to):
      er = 6366.707

      latFrom = deg2rad(self.lat)
      latTo   = deg2rad(to.lat)
      lngFrom = deg2rad(self.lng)
      lngTo   = deg2rad(to.lng)

      x1 = er * cos(lngFrom) * sin(latFrom)
      y1 = er * sin(lngFrom) * sin(latFrom)
      z1 = er * cos(latFrom)

      x2 = er * cos(lngTo) * sin(latTo)
      y2 = er * sin(lngTo) * sin(latTo)
      z2 = er * cos(latTo)

      d = acos(sin(latFrom)*sin(latTo) + cos(latFrom)*cos(latTo)*cos(lngTo-lngFrom)) * er
      
      return d

    
    """
     * Convert this LatLng object from OSGB36 datum to WGS84 datum.
    """
    def OSGB36ToWGS84(self):
        airy1830 = RefEll(6377563.396, 6356256.909)
        a        = airy1830.maj
        b        = airy1830.min
        eSquared = airy1830.ecc
        phi = deg2rad(self.lat)
        lam = deg2rad(self.lng)
        v = a / (sqrt(1 - eSquared * sinSquared(phi)))
        H = 0 # height
        x = (v + H) * cos(phi) * cos(lam)
        y = (v + H) * cos(phi) * sin(lam)
        z = ((1 - eSquared) * v + H) * sin(phi)

        tx =        446.448
        ty =       -124.157
        tz =        542.060
        s  =         -0.0000204894
        rx = deg2rad( 0.00004172222)
        ry = deg2rad( 0.00006861111)
        rz = deg2rad( 0.00023391666)

        xB = tx + (x * (1 + s)) + (-rx * y)     + (ry * z)
        yB = ty + (rz * x)      + (y * (1 + s)) + (-rx * z)
        zB = tz + (-ry * x)     + (rx * y)      + (z * (1 + s))

        wgs84 = RefEll(6378137.000, 6356752.3141)
        a        = wgs84.maj
        b        = wgs84.min
        eSquared = wgs84.ecc

        lambdaB = rad2deg(atan(yB / xB))
        p = sqrt((xB * xB) + (yB * yB))
        phiN = atan(zB / (p * (1 - eSquared)))

        for i in range(1,10):
            v = a / (sqrt(1 - eSquared * sinSquared(phiN)))
            phiN1 = atan((zB + (eSquared * v * sin(phiN))) / p)
            phiN = phiN1

        phiB = rad2deg(phiN)

        self.lat = phiB
        self.lng = lambdaB
    
    
    """
     * Convert this LatLng object from WGS84 datum to OSGB36 datum.
    """
    def WGS84ToOSGB36(self):
        wgs84 = RefEll(6378137.000, 6356752.3141)
        a        = wgs84.maj
        b        = wgs84.min
        eSquared = wgs84.ecc
        phi = deg2rad(self.lat)
        lam = deg2rad(self.lng)
        v = a / (sqrt(1 - eSquared * sinSquared(phi)))
        H = 0; # height
        x = (v + H) * cos(phi) * cos(lam)
        y = (v + H) * cos(phi) * sin(lam)
        z = ((1 - eSquared) * v + H) * sin(phi)

        tx =       -446.448
        ty =        124.157
        tz =       -542.060
        s  =          0.0000204894
        rx = deg2rad(-0.00004172222)
        ry = deg2rad(-0.00006861111)
        rz = deg2rad(-0.00023391666)

        xB = tx + (x * (1 + s)) + (-rx * y)     + (ry * z)
        yB = ty + (rz * x)      + (y * (1 + s)) + (-rx * z)
        zB = tz + (-ry * x)     + (rx * y)      + (z * (1 + s))

        airy1830 = RefEll(6377563.396, 6356256.909)
        a        = airy1830.maj
        b        = airy1830.min
        eSquared = airy1830.ecc

        lambdaB = rad2deg(atan(yB / xB))
        p = sqrt((xB * xB) + (yB * yB))
        phiN = atan(zB / (p * (1 - eSquared)))
      
        for i in range(1,10):
            v = a / (sqrt(1 - eSquared * sinSquared(phiN)))
            phiN1 = atan((zB + (eSquared * v * sin(phiN))) / p)
            phiN = phiN1
 
        phiB = rad2deg(phiN)

        self.lat = phiB
        self.lng = lambdaB
    
    
    """
     * Convert this LatLng object into an OSGB grid reference. Note that this
     * function does not take into account the bounds of the OSGB grid -
     * beyond the bounds of the OSGB grid, the resulting OSRef object has no
     * meaning
     *
     * @return the converted OSGB grid reference
    """
    def toOSRef(self):
        airy1830 = RefEll(6377563.396, 6356256.909)
        OSGB_F0  = 0.9996012717
        N0       = -100000.0
        E0       = 400000.0
        phi0     = deg2rad(49.0)
        lambda0  = deg2rad(-2.0)
        a        = airy1830.maj
        b        = airy1830.min
        eSquared = airy1830.ecc
        phi = deg2rad(self.lat)
        lam = deg2rad(self.lng)
        E = 0.0
        N = 0.0
        n = (a - b) / (a + b)
        v = a * OSGB_F0 * pow(1.0 - eSquared * sinSquared(phi), -0.5)
        rho = a * OSGB_F0 * (1.0 - eSquared) * pow(1.0 - eSquared * sinSquared(phi), -1.5)
        etaSquared = (v / rho) - 1.0
        M = \
            (b * OSGB_F0) \
              * (((1 + n + ((5.0 / 4.0) * n * n) + ((5.0 / 4.0) * n * n * n)) \
                * (phi - phi0)) \
                - (((3 * n) + (3 * n * n) + ((21.0 / 8.0) * n * n * n)) \
                  * sin(phi - phi0) \
                  * cos(phi + phi0)) \
                + ((((15.0 / 8.0) * n * n) + ((15.0 / 8.0) * n * n * n)) \
                  * sin(2.0 * (phi - phi0)) \
                  * cos(2.0 * (phi + phi0))) \
                - (((35.0 / 24.0) * n * n * n) \
                  * sin(3.0 * (phi - phi0)) \
                  * cos(3.0 * (phi + phi0))))
        I = M + N0
        II = (v / 2.0) * sin(phi) * cos(phi)
        III = \
            (v / 24.0) \
            * sin(phi) \
            * pow(cos(phi), 3.0) \
            * (5.0 - tanSquared(phi) + (9.0 * etaSquared))
            
        IIIA = \
            (v / 720.0) \
            * sin(phi) \
            * pow(cos(phi), 5.0) \
            * (61.0 - (58.0 * tanSquared(phi)) + pow(tan(phi), 4.0))
            
        IV = v * cos(phi)
        V = (v / 6.0) * pow(cos(phi), 3.0) * ((v / rho) - tanSquared(phi))
        VI = \
            (v / 120.0) \
            * pow(cos(phi), 5.0) \
            * (5.0 \
                - (18.0 * tanSquared(phi)) \
                + (pow(tan(phi), 4.0)) \
                + (14 * etaSquared) \
                - (58 * tanSquared(phi) * etaSquared))

        N = \
            I \
            + (II * pow(lam - lambda0, 2.0)) \
            + (III * pow(lam - lambda0, 4.0)) \
            + (IIIA * pow(lam - lambda0, 6.0))
        E = \
            E0 \
            + (IV * (lam - lambda0)) \
            + (V * pow(lam - lambda0, 3.0)) \
            + (VI * pow(lam - lambda0, 5.0))

        return OSRef(E, N)
    
    
    """
     * Convert a latitude and longitude to an UTM reference
     *
     * @return the converted UTM reference
    """
    def toUTMRef(self):
        wgs84 = RefEll(6378137, 6356752.314)
        UTM_F0   = 0.9996
        a = wgs84.maj
        eSquared = wgs84.ecc
        longitude = self.lng
        latitude = self.lat

        latitudeRad = latitude * (pi() / 180.0)
        longitudeRad = longitude * (pi() / 180.0)
        longitudeZone = (int) ((longitude + 180.0) / 6.0) + 1

        # Special zone for Norway
        if (latitude >= 56.0 \
                and latitude < 64.0 \
                and longitude >= 3.0 \
                and longitude < 12.0):
            longitudeZone = 32

        # Special zones for Svalbard
        if latitude >= 72.0 and latitude < 84.0:
            if longitude >= 0.0 and longitude < 9.0:
                longitudeZone = 31
            elif longitude >= 9.0 and longitude < 21.0:
                longitudeZone = 33
            elif longitude >= 21.0 and longitude < 33.0:
                longitudeZone = 35
            elif longitude >= 33.0 and longitude < 42.0:
                longitudeZone = 37
      
        longitudeOrigin = (longitudeZone - 1) * 6 - 180 + 3
        longitudeOriginRad = longitudeOrigin * (pi() / 180.0)

        UTMZone = getUTMLatitudeZoneLetter(latitude)

        ePrimeSquared = (eSquared) / (1 - eSquared)

        n = a / sqrt(1 - eSquared * sin(latitudeRad) * sin(latitudeRad))
        t = tan(latitudeRad) * tan(latitudeRad)
        c = ePrimeSquared * cos(latitudeRad) * cos(latitudeRad)
        A = cos(latitudeRad) * (longitudeRad - longitudeOriginRad)

        M = \
            a \
            * ((1 \
            - eSquared / 4 \
            - 3 * eSquared * eSquared / 64 \
            - 5 * eSquared * eSquared * eSquared / 256) \
            * latitudeRad \
            - (3 * eSquared / 8 \
              + 3 * eSquared * eSquared / 32 \
              + 45 * eSquared * eSquared * eSquared / 1024) \
              * sin(2 * latitudeRad) \
            + (15 * eSquared * eSquared / 256 \
              + 45 * eSquared * eSquared * eSquared / 1024) \
              * sin(4 * latitudeRad) \
            - (35 * eSquared * eSquared * eSquared / 3072) \
              * sin(6 * latitudeRad))

        UTMEasting = \
            float(UTM_F0 \
            * n \
            * (A \
                + (1 - t + c) * pow(A, 3.0) / 6 \
                + (5 - 18 * t + t * t + 72 * c - 58 * ePrimeSquared) \
                * pow(A, 5.0) \
                / 120) \
            + 500000.0) \

        UTMNorthing = \
            float(UTM_F0 \
                * (M \
                + n \
                * tan(latitudeRad) \
                * (A * A / 2 \
                    + (5 - t + (9 * c) + (4 * c * c)) * pow(A, 4.0) / 24 \
                    + (61 - (58 * t) + (t * t) + (600 * c) - (330 * ePrimeSquared)) \
                    * pow(A, 6.0) \
                    / 720)))

        # Adjust for the southern hemisphere
        if (latitude < 0):
            UTMNorthing += 10000000.0

        return UTMRef(UTMEasting, UTMNorthing, UTMZone, longitudeZone)


# =================================================================== OSRef

# References given with OSRef are accurate to 1m.
class OSRef:

    easting = 0
    northing = 0


    """
     * Create a new OSRef object representing an OSGB grid reference. Note
     * that the parameters for this constructor require eastings and
     * northings with 1m accuracy and need to be absolute with respect to
     * the whole of the British Grid. For example, to create an OSRef
     * object from the six-figure grid reference TG514131, the easting would
     * be 651400 and the northing would be 313100.
     * 
     * Grid references with accuracy greater than 1m can be represented
     * using floating point values for the easting and northing. For example,
     * a value representing an easting or northing accurate to 1mm would be
     * given as 651400.0001.
     *
     * @param easting the easting of the reference (with 1m accuracy)
     * @param northing the northing of the reference (with 1m accuracy)
    """
    def __init__(self,easting, northing):
      self.easting  = easting
      self.northing = northing


    """
     * Convert this grid reference into a string showing the exact values
     * of the easting and northing.
     *
     * @return
    """
    def __str__(self):
        return "(%f,%f)" % (self.easting,self.northing)


    """
    * Take a string formatted as a six-figure OS grid reference (e.g.
    * "TG514131") and return a reference to an OSRef object that represents
    * that grid reference. The first character must be H, N, S, O or T.
    * The second character can be any uppercase character from A through Z
    * excluding I.
    *
    * @param ref
    * @return
    * @since 2.1
    """
    @staticmethod
    def fromSixFigureString(ref):
        char1 = ref[0:1]
        char2 = ref[1:2]
        east  = int(ref[2:5]) * 100
        north = int(ref[5:8]) * 100
        if (char1 == 'H'):
          north += 1000000
        elif (char1 == 'N'):
          north += 500000
        elif (char1 == 'O'):
          north += 500000
          east  += 500000
        elif (char1 == 'T'):
          east += 500000
        
        char2ord = ord(char2)
        if (char2ord > 73): char2ord-=1; # Adjust for no I
        nx = ((char2ord - 65) % 5) * 100000
        ny = (4 - floor((char2ord - 65) / 5)) * 100000
        return OSRef(east + nx, north + ny)
        
    """
     * Convert this grid reference into a string using a standard six-figure
     * grid reference including the two-character designation for the 100km
     * square. e.g. TG514131. 
     *
     * @return
    """
    def toSixFigureString(self):
        hundredkmE = floor(self.easting / 100000)
        hundredkmN = floor(self.northing / 100000)
        firstLetter = ""
        if (hundredkmN < 5):
            if (hundredkmE < 5):
                firstLetter = "S"
            else:
                firstLetter = "T"
            
        elif (hundredkmN < 10):
            if (hundredkmE < 5):
                firstLetter = "N"
            else:
                firstLetter = "O"
        else:
            firstLetter = "H"

        secondLetter = ""
        index = int(65 + ((4 - (hundredkmN % 5)) * 5) + (hundredkmE % 5))
        ti = index
        if (index >= 73): index+=1
        secondLetter = chr(index)

        e = round((self.easting - (100000 * hundredkmE)) / 100)
        n = round((self.northing - (100000 * hundredkmN)) / 100)

        return "%s%s%03d%03d" % (firstLetter, secondLetter, e, n)


    """
     * Convert this grid reference into a latitude and longitude
     *
     * @return
    """
    def toLatLng(self):
        airy1830 = RefEll(6377563.396, 6356256.909)
        OSGB_F0  = 0.9996012717
        N0       = -100000.0
        E0       = 400000.0
        phi0     = deg2rad(49.0)
        lambda0  = deg2rad(-2.0)
        a        = airy1830.maj
        b        = airy1830.min
        eSquared = airy1830.ecc
        phi      = 0.0
        lam   = 0.0
        E        = self.easting
        N        = self.northing
        n        = (a - b) / (a + b)
        M        = 0.0
        phiPrime = ((N - N0) / (a * OSGB_F0)) + phi0
      
        while True:
            M = \
                (b * OSGB_F0) \
                    * (((1 + n + ((5.0 / 4.0) * n * n) + ((5.0 / 4.0) * n * n * n)) \
                    * (phiPrime - phi0)) \
                    - (((3 * n) + (3 * n * n) + ((21.0 / 8.0) * n * n * n)) \
                        * sin(phiPrime - phi0) \
                        * cos(phiPrime + phi0)) \
                    + ((((15.0 / 8.0) * n * n) + ((15.0 / 8.0) * n * n * n)) \
                        * sin(2.0 * (phiPrime - phi0)) \
                        * cos(2.0 * (phiPrime + phi0))) \
                    - (((35.0 / 24.0) * n * n * n) \
                        * sin(3.0 * (phiPrime - phi0)) \
                        * cos(3.0 * (phiPrime + phi0)))) \
                        
            phiPrime += (N - N0 - M) / (a * OSGB_F0)
            
            if ((N - N0 - M) >= 0.001):
                break
                
        v = a * OSGB_F0 * pow(1.0 - eSquared * sinSquared(phiPrime), -0.5)
        rho = \
            a \
            * OSGB_F0 \
            * (1.0 - eSquared) \
            * pow(1.0 - eSquared * sinSquared(phiPrime), -1.5)
            
        etaSquared = (v / rho) - 1.0
        VII = tan(phiPrime) / (2 * rho * v)
        VIII = \
            (tan(phiPrime) / (24.0 * rho * pow(v, 3.0))) \
            * (5.0 \
                + (3.0 * tanSquared(phiPrime)) \
                + etaSquared \
                - (9.0 * tanSquared(phiPrime) * etaSquared))
        IX = \
            (tan(phiPrime) / (720.0 * rho * pow(v, 5.0))) \
            * (61.0 \
                + (90.0 * tanSquared(phiPrime)) \
                + (45.0 * tanSquared(phiPrime) * tanSquared(phiPrime)))
                
        X = sec(phiPrime) / v
        XI = \
            (sec(phiPrime) / (6.0 * v * v * v)) \
            * ((v / rho) + (2 * tanSquared(phiPrime)))
        XII = \
            (sec(phiPrime) / (120.0 * pow(v, 5.0))) \
            * (5.0 \
                + (28.0 * tanSquared(phiPrime)) \
                + (24.0 * tanSquared(phiPrime) * tanSquared(phiPrime)))
        XIIA = \
            (sec(phiPrime) / (5040.0 * pow(v, 7.0))) \
            * (61.0 \
                + (662.0 * tanSquared(phiPrime)) \
                + (1320.0 * tanSquared(phiPrime) * tanSquared(phiPrime)) \
                + (720.0 \
                * tanSquared(phiPrime) \
                * tanSquared(phiPrime) \
                * tanSquared(phiPrime)))
        phi = \
            phiPrime \
            - (VII * pow(E - E0, 2.0)) \
            + (VIII * pow(E - E0, 4.0)) \
            - (IX * pow(E - E0, 6.0))
        lam = \
            lambda0 \
            + (X * (E - E0)) \
            - (XI * pow(E - E0, 3.0)) \
            + (XII * pow(E - E0, 5.0)) \
            - (XIIA * pow(E - E0, 7.0))
 
        return LatLng(rad2deg(phi), rad2deg(lam))

        
  # ================================================================== UTMRef

class UTMRef:

    easting = 0
    northing = 0
    latZone = 0
    lngZone = 0


    """
     * Create a new object representing a UTM reference.
     *
     * @param easting
     * @param northing
     * @param latZone
     * @param lngZone
    """
    def __init__(self,easting, northing, latZone, lngZone):
      self.easting  = easting
      self.northing = northing
      self.latZone  = latZone
      self.lngZone  = lngZone


    """
     * Return a string representation of this UTM reference
     *
     * @return
    """
    def toString(self):
      return "%s%s %s%s" % map(str,(self.lngZone, self.latZone, self.easting, self.northing))
    
    """
     * Convert this UTM reference to a latitude and longitude
     *
     * @return the converted latitude and longitude
    """
    def toLatLng(self):
      wgs84 = RefEll(6378137, 6356752.314)
      UTM_F0   = 0.9996
      a = wgs84.maj
      eSquared = wgs84.ecc
      ePrimeSquared = eSquared / (1.0 - eSquared)
      e1 = (1 - sqrt(1 - eSquared)) / (1 + sqrt(1 - eSquared))
      x = self.easting - 500000.0;
      y = self.northing
      zoneNumber = self.lngZone
      zoneLetter = self.latZone

      longitudeOrigin = (zoneNumber - 1.0) * 6.0 - 180.0 + 3.0

      # Correct y for southern hemisphere
      if ((ord(zoneLetter) - ord("N")) < 0):
        y -= 10000000.0

      m = y / UTM_F0
      mu = \
        m \
          / (a \
            * (1.0 \
              - eSquared / 4.0 \
              - 3.0 * eSquared * eSquared / 64.0 \
              - 5.0 \
                * pow(eSquared, 3.0) \
                / 256.0))

      phi1Rad = \
        mu \
          + (3.0 * e1 / 2.0 - 27.0 * pow(e1, 3.0) / 32.0) * sin(2.0 * mu) \
          + (21.0 * e1 * e1 / 16.0 - 55.0 * pow(e1, 4.0) / 32.0) \
            * sin(4.0 * mu) \
          + (151.0 * pow(e1, 3.0) / 96.0) * sin(6.0 * mu)

      n = \
        a \
          / sqrt(1.0 - eSquared * sin(phi1Rad) * sin(phi1Rad))
      t = tan(phi1Rad) * tan(phi1Rad)
      c = ePrimeSquared * cos(phi1Rad) * cos(phi1Rad)
      r = \
        a \
          * (1.0 - eSquared) \
          / pow( \
            1.0 - eSquared * sin(phi1Rad) * sin(phi1Rad), \
            1.5)
      d = x / (n * UTM_F0)

      latitude = (
        phi1Rad
          - (n * tan(phi1Rad) / r)
            * (d * d / 2.0
              - (5.0
                + (3.0 * t)
                + (10.0 * c)
                - (4.0 * c * c)
                - (9.0 * ePrimeSquared))
                * pow(d, 4.0)
                / 24.0
              + (61.0
                + (90.0 * t)
                + (298.0 * c)
                + (45.0 * t * t)
                - (252.0 * ePrimeSquared)
                - (3.0 * c * c))
                * pow(d, 6.0)
                / 720.0)) * (180.0 / pi())

      longitude = longitudeOrigin + (
        (d
          - (1.0 + 2.0 * t + c) * pow(d, 3.0) / 6.0
          + (5.0
            - (2.0 * c)
            + (28.0 * t)
            - (3.0 * c * c)
            + (8.0 * ePrimeSquared)
            + (24.0 * t * t))
            * pow(d, 5.0)
            / 120.0)
          / cos(phi1Rad)) * (180.0 / pi())

      return LatLng(latitude, longitude)


  # ================================================================== RefEll

class RefEll:

    maj = 0
    min = 0
    ecc = 0


    """
     * Create a new RefEll object to represent a reference ellipsoid
     *
     * @param maj the major axis
     * @param min the minor axis
    """
    def __init__(self, maj, min):
      self.maj = maj
      self.min = min
      self.ecc = ((maj * maj) - (min * min)) / (maj * maj)


  # ================================================== Mathematical Functions

def sinSquared(x):
    return sin(x) * sin(x)

def cosSquared(x):
    return cos(x) * cos(x)

def tanSquared(x):
    return tan(x) * tan(x)

def sec(x):
    return 1.0 / cos(x)
  
_d2r = atan(1.0) * 4 / 180.0
  
def deg2rad(x): return x * _d2r
def rad2deg(x): return x / _d2r
  
"""
   *  Work out the UTM latitude zone from the latitude
   *
   * @param latitude
   * @return
"""
def getUTMLatitudeZoneLetter(latitude):
    if ((84 >= latitude) and (latitude >= 72)): return "X"
    elif (( 72 > latitude) and (latitude >=  64)): return "W"
    elif (( 64 > latitude) and (latitude >=  56)): return "V"
    elif (( 56 > latitude) and (latitude >=  48)): return "U"
    elif (( 48 > latitude) and (latitude >=  40)): return "T"
    elif (( 40 > latitude) and (latitude >=  32)): return "S"
    elif (( 32 > latitude) and (latitude >=  24)): return "R"
    elif (( 24 > latitude) and (latitude >=  16)): return "Q"
    elif (( 16 > latitude) and (latitude >=   8)): return "P"
    elif ((  8 > latitude) and (latitude >=   0)): return "N"
    elif ((  0 > latitude) and (latitude >=  -8)): return "M"
    elif (( -8 > latitude) and (latitude >= -16)): return "L"
    elif ((-16 > latitude) and (latitude >= -24)): return "K"
    elif ((-24 > latitude) and (latitude >= -32)): return "J"
    elif ((-32 > latitude) and (latitude >= -40)): return "H"
    elif ((-40 > latitude) and (latitude >= -48)): return "G"
    elif ((-48 > latitude) and (latitude >= -56)): return "F"
    elif ((-56 > latitude) and (latitude >= -64)): return "E"
    elif ((-64 > latitude) and (latitude >= -72)): return "D"
    elif ((-72 > latitude) and (latitude >= -80)): return "C"
    else: return 'Z'
